Asymptotics of Reaction-Diffusion Fronts with One 
Static and One Diffusing Reactant 



Martin Z. Bazant 

Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139 

H. A. Stone 

Division of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138 

June 8, 2000 



Abstract 

The long-time behavior of a reaction-diffusion front between one static {e.g. porous 
sohd) reactant A and one initiaUy separated diffusing reactant B is analyzed for the 
mean-field reaction-rate density R{pa,Pb) = kp^p%. A uniformly valid asymptotic 
approximation is constructed from matched self-similar solutions in a "reaction front" 
(of width tt; ~ where R ^ t^ enters the dominant balance) and a "diffusion layer" 
(of width W ~ where R is negligible). The limiting solution exists if and only if 
m, n > 1, in which case the scaling exponents are uniquely given by a = {m — \)/2{m+l) 
and (5 = m/{m + 1). In the diffusion layer, the common ad hoc approximation of 
neglecting reactions is given mathematical justification, and the exact transient decay of 
the reaction rate is derived. The physical effects of higher-order kinetics (m, n > 1), such 
as the broadening of the reaction front and the slowing of transients, are also discussed. 
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1 Introduction 



In the decade that has passed since the pioneering analytical study of Galfi and Racz [|ll, there 
has emerged a substantial body of research devoted to experimental 0, ^, computational 0, 
H) S 0' B B HOI analytical [0, [l^, |T3|, |14|, |l^, |T6| studies of reaction-diffusion systems 
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with two initially separated, diffusing species A and B reacting to produce an inert product 
C according to the chemical formula 



m'A + n'B ^ C (inert), (1) 

where m' and n', the stoichiometric coefficients, are positive integers. Theoretical studies have 
focused almost exclusively on the "one- dimensional" case of an infinite, fiat reaction front 
between two regions of homogeneous composition of either A or B (see Fig. |I|). This idealized 
situation is believed to capture much of the essential physics of reaction fronts commonly 



observed in various chemical |T^, |T8[ and biological |T9|, ^ systems. 

The standard continuum model for such a one-dimensional reaction front involves a pair 
of nonhnear partial differential equations ^, |TT], [T^ |T^ 



^ = D^^-m'R{p^,ps) (2a) 

^ = D,^-n'Rip,,p,), (2b) 

subject to the boundary conditions 

Pa(-oo,T) = 0, pAioo,T) = p% pb{-oo,T) = p%, pb(oo,T) = (3) 

and the initial conditions 

Pa{X, 0) = p\H{X), pb{X, 0) = p°sH{-X) (4) 

where p^(X, T) and pb{X,T) are the concentrations, and Db the diffusion coefficients of 
A and B, respectively, > and p^ > are constants, H{X) is the Heaviside unit step func- 
tion and R{pa, Pb) is the reaction rate density for production of species C. (Note that upper- 
case letters denote quantities with dimensions, e.g. X and T for space and time, respectively. 
Lower-case letters for the corresponding dimensionless quantities are introduced in section |2.1| .) 
The reactants are completely separated at first according to (||), but for T > they diffuse 
together and react, which decreases the concentrations wherever pa{X,T)pb{X,T) > 0. Dif- 
fusion acts to replenish any depleted regions. As a result the system develops a localized, 
moving region, the "reaction front," where the reaction rate R{pa,Pb) is greatest and which 
is fed by diffusion from the distant particle reservoirs described by the boundary conditions. 
The dynamics of this reaction front are described by the long-time asymptotics of the nonlinear 
initial-boundary- value problem (2)-(^. 

The nonlinear reaction term R{pa, P_b) is usually assumed to have the form of a power law 

R{pA,PB)=kp'^pl,, (5) 

where /c is a rate constant, and m and n are respectively the "kinetic orders" of A and B in the 
reaction |^. For a one-step reaction with sufficient mixing (see below) m = m' and n = n', 
but for more complex, multi-step reactions m and n are determined by the stoichiometric 
coefficients of the (often unknown) rate-limiting step. Although m and n are usually taken 
to be positive integers, non- integer values of m and n can arise in certain situations |]18|. We 
will see that a well-defined reaction front exists for any real numbers m, n > 1, but not for 
m < 1 or n < 1. 
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Technically, by assuming in (2) that the reaction rate R depends only on the average 
local concentrations (and not on any fluctuations or many-body effects) we have made the 



"mean-field approximation" |21|. In low- dimensional systems, such as ion channels [d = 1) or 
catalytic surfaces (rf = 2), the mean-field approximation can break down because the reacting 
particles cannot mix efficiently enough, but as the dimension of the system is increased above 
a certain "upper critical dimension" dc, such statistical anomalies disappear. For two diffusing 
reactants with a simple one-step reaction it is known [^, |^, ^ that dc = 2/{m + n — l). Since 
dc < 2 for 171,11 > 1 the mean-field approximation should be perfectly valid in the usual case 
d = 3, which is consistent with experimental findings 0. 

In contrast to the case of two diffusing reactants described above, relatively little is 



known ^, |T^, [T^, 16 1 about the case of one diffusing reactant {Da > 0) and one static 
reactant {Db = 0)- This situation, depicted schematically in Fig. |I], describes the corrosion of 
a porous solid B saturated with a fluid solvent and exposed to an initially separated colloidal 
reactant A, as shown in recent electrochemical experiments (described below) 0. Jiang and 
Ebner first pointed out (for m = n = 1) that setting = in (2) is a non-trivial, i.e. 
singular, limit leading to different long-time behavior than in the case of Db > (no matter 
how small), which they explained with simple scaling arguments supported by Monte Carlo 
computer simulations. For an analytical description of such one-dimensional diffusion with 
one static reactant, we adopt the power-law form of the reaction term and study the coupled 
equations 

~&f ^ ^A^-rnkp^^pl (6a) 
dpB 



dT 



n'kp^pl- (6b) 



In the simplest case m = n = 1, the initial-boundary- value problem (^-(6) has been solved 
numerically by Havlin et al. ||^ and analyzed in the limit of "long times" T — ^> oo by Koza ||14|| , 
using various asymptotic approximations introduced by Galfi and Racz . Rigorous analysis 
has been reported in the analogous limit of "fast reactions" — > oo by Hilhorst et al. , but 



these authors only address the behavior at the diffusive length scale X oc vT (see section 3.4 



below) and do not consider the structure of the reaction front studied by Koza [|14|, which 



is of primary interest here. Hilhorst et al. have also recently considered the effect of a more 



general reaction term at the diffusive scale [If:], but the present work appears to be the first 



to analyze the nontrivial effect of changing reaction orders at the reactive length scale (see 
below) in the general case (m, n > 1) with one static reactant. 

The relevance of (6) for a given porous-solid corrosion system rests on several key assump- 
tions that are less obviously satisfied a priori than in the case of two diffusing reactants. First, 
the solid matrix containing the static reactant B must be sufficiently porous that the mov- 
ing reactant A can diffuse freely to the exposed surfaces with an effective diffusion constant 
(averaged over many pores) comparable to that in the bulk solvent. The concentration of A 
must also be dilute enough that D^ is constant. Another reason that the concentrations of 
A and B must be dilute is that the inert product C must be created in small enough quan- 
tities that its presence does not affect the reaction dynamics (e.g. by inhibiting diffusion or 
initiating convection). Finally, one might worry about the breakdown of the mean- field ap- 
proximation since the (possibly fractal) pore structure may influence statistical averaging. For 
example, it is known that fluctuations alter the reaction-front dynamics when the diffusion is 
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confined to a percolating cluster in two dimensions In spite of these concerns, however, 
the one-dimensional mean-field model (6) can in fact describe certain corrosion systems. 

An important motivation for the present analytical study is afforded by the recent exper- 
iments of Leger et al. [^], which are the first to examine in detail the case of one static and 
one diffusing reactant. These experiments involve the corrosion of ramified copper electrode- 
posits exposed to a cupric chloride electrolyte to produce cuprous chloride crystallites via the 
reaction 

CuCla (aq) + Cu (solid) ^ 2 CuCl (solid) (7) 

immediately following electrodeposition. It is found that the long-time behavior of (6) with 
m = 1 matches the experimentally observed front speed and concentration profile of diffusing 
reactant (CUCI2) rather well, in spite of the complex fractal geometry of the electrodeposits 
and the presence of the inert product (CuCl) [||. Since the reaction rate and the concentration 
of the static species (Cu) are not directly measured, however, the interpretation of these kinds 
of corrosion experiments can be aided by the analysis presented here of the mean-field model 
with m,n > 1. 

There is an extensive mathematical literature [p!7| , |19|, |23|, on the subject of single 
reaction-diffusion equations of the general form 

which arise in many applications (e.g. chemical reactions, combustion and population dynam- 
ics). A common theme in these studies is the appearance of two distinct (time-dependent) 
length scales in the intermediate asymptotic regime {t — > 00) which correspond to either 
"weakly nonlinear behavior", where it has been established in many cases that the reaction 
term is negligible and the dynamics are purely diffusive, or "strongly nonlinear behavior", 
where the reaction and diffusion terms balance (in the nomenclature of Gmira and Veron ||23|| ). 
This separation of scales also arises in coupled systems of reaction-diffusion equations like (2), 
but owing to their greater complexity, much less rigorous analysis has been reported. In the 
case of two diffusing reactants, Galfi and Racz [|T| pointed out that if the diffusion constants 
are the same. Da = Db, then the difference in concentrations pA — Pb obeys a pure diffusion 
equation which can be easily integrated, thereby reducing the coupled system (2) to a single 
equation with the form of (P). Another simplification occurs if also = p% in which case 
the reaction front is perfectly symmetric and does not move. In this simplified case with 
m = n = 1, Schenkel et al. were able to prove that the asymptotic solution of Galfi 



and Racz 0], which combines different approximations at the diffusive and reactive scales, is 
approached uniformly as T — > cxd starting from the initial conditions of (^, and they also 
reported rigorous bounds on the transient decay to the asymptotic solution. Recently, van 
Baalen et al. have extended this analysis to the case of symmetric, high-order reactions 
m = n > 3, where the reaction-front scaling is altered. 

The analyses of Refs. |Tl|-|T^ represent an important contribution because, at least in the 



case Da = Db, Pa = P°b m = n, they provide a rigorous mathematical justification for 
various ad hoc assumptions introduced by Galfi and Racz ^ |TB|] to describe the local structure 
of the reaction front which have otherwise been validated only by numerical simulations. 
Unfortunately, however, since the analysis in Refs. |ll|, |12| relies on a comparison principle 



for single parabolic equations |2^, ^ of the form (§), it does not (as the authors indicate) 
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appear to be applicable when Db (which also leads to a moving reaction front). Van 

Baalen et al. [|1^] also remark that their analysis is not easily extended to certain intermediate 



reaction orders (1 < m = < 3). These difficulties are reflected in Koza's recent studies of the 



general cases Da > Db > [0 and Db = fT^, in which several ad hoc (but reasonable) 
approximations are made and transients are ignored. 

In the present article, the long-time asymptotics of the initial-boundary- value problem (|^)- 
(6) are studied. This special case of (2) is more tractable analytically than the general case 
because (6b) can be integrated exactly in time, thereby reducing the coupled system to a single 
integro-partial differential equation. This useful simplification is presented in section ^ where 
the problem is recast in a dimensionless form. It is also noted that similarity solutions are 
expected to exist because there is no natural length or time scale in the problem [24, 26[ . 
Although it may be possible to prove that the system actually approaches such a self-similar 
solution starting from the prescribed initial conditions, we instead pursue the more modest 
goal of proving that if an asymptotic similarity solution exists, it must have a certain unique 
form, i.e. we explore the consequences of the "quasi-stationary approximation" |T3|, In 



section ^, the similarity solution is systematically derived, and it is shown that a "diffusion 
layer" (where the reaction term is dominated by the diffusion term) with different scaling 
properties than the "reaction front" (where the reaction and diffusion terms balance) must 
exist to satisfy the boundary conditions. In section ^, the transient decay of the reaction 
rate in the diffusion layer is analyzed, thereby proving a posteriori that the reaction term 
can indeed be neglected in the dominant balance. In section ^, a uniformly- valid, asymptotic 
approximation is constructed by matching the self-similar forms in the two different regions. 
Finally, in section ^ some general physical conclusions are drawn from the analysis, and in the 
Epilogue certain similarities are discussed between this work and the literature on combustion 
waves. (Note that section | is more technical and may be skipped in a first reading.) 

2 Preliminaries 

2.1 Dimensionless Formulation 

With the definitions, 

t^m'k{p'Xr~\p''BrT, x^X^m'k{p\Y^-\p%Y/DA (9a) 
a(x, t) = pa{X, T)/p% b{x, t) = pb{X, T)/p%, (9b) 

the initial-boundary-value problem (|^)-(6) may be expressed in a dimensionless form 

db 

— = -qaJ^b"^. (10b) 

a(cx3,t) = 1, b{oo,t) = 0, a{—oo,t) = 0, 6(— oo,t) = 1, (10c) 
a(x, 0) = H{x), b{x, 0) = H{-x) (lOd) 

which involves only one dimensionless parameter: 

q^^. (11) 
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Note that the dimensionless problem (|T^) depends only upon the initial concentrations p\ 
and rho% and the stochiometric coefficients m! and n' through the parameter g; the reaction 
rate k are the diffusion constant Da simply set the natural scales for length and time. From 
(|9aD, we see that the limit of "fast reactions" k oo (with X and T fixed) corrseponds to 
the limit of long (dimensionless) times t ^ oo at the diffusive scale x oc \/t. (See Ref. [|15| for 
another discussion of this correspondence of limits.) 



2.2 The Governing Integro- Partial Differential Equation 



The statement of the problem (|TD|) will be used in deriving the asymptotic similarity solution 
below, but for the transient analysis described in section ^ it will be convenient to ffist 
integrate (10b) exactly in time. Note that (10b) and ( |10d| ) imply that b{x,t) = for x > at 
all times t > 0, which reflects the fact that species B cannot diffuse out of its initial region. 
For X < 0, we integrate (10b) using the initial condition (|10d|) to express b{x,t) as 



b{x,t) 



'q<j)m(x,t) 



^ [l + g(n-l)0„(x,t)]-^/("-^) 
which involves the time-integral of a(x, t)™': 



if 72 = 1 



a; < 



(12) 



a x, r 



'dr. 



(13) 



(Note that a partial differential equation satisfied by 0i(x, t) is given in Ref. JlSl.) Substituting 
for b{x, t) in ( |10a| ), we obtain a single, nonlinear integro-partial differential equation for a(x, t), 
either 

= ^!^^^ _ i7(x)a(x,t)-e-''/o 



if n = 1 or 



da{x,t) d'^a{x,t) 



H{x)a{x, ty 



dt 



1 + q{n — 1) /q a(x, T)"^dT 



n/(n— 1) 



(14) 
(15) 



if n ^ 1. 

Although these equations involve only one unknown function a{x,t), they are somewhat 
unwieldy, so we will ffist seek long-time {t oo) asymptotic solutions to the coupled system 
( PUj) in section ^. The time-dependent properties of (H) and (^) will be studied in section]^. 
Before proceeding, however, we digress to show that physically meaningful solutions exist only 
if n > 1. Later in the analysis we will see that m > 1 is required as well. 



2.3 A Reaction Front Does Not Exist if n < 1 

Consider any point Xo < 0. Since species A diffuses to Xq from a reservoir of constant concen- 
tration (a(oo,t) = 1) while species B is removed by reactions without ever being replenished 
{db{xo,t)/dt < for all t > 0), it is clear that after long times a{xo,t) must eventually differ 
from zero. Therefore, there exists some a*(xo) > and to > such that a{xo,t) > a*(xo) for 
all t > to- This implies (pmixo.t) > a*(xo)™ ■ (t — to) from (p^) and thus b{xo,t) —>■ from 
(|T2p since q > 0, but a singularity arises if < 1: The concentration of static reactant b{xo,t) 
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vanishes at some finite time ti given by g(l — n)0m(xo,ti) = 1 (which exists because 4>{xo,t) 
is continuous, (j)mixo,0) = and 0^(^^0,00) = 00). For t > ti, Eq. (|1^) predicts imaginary, 
negative, or diverging solutions for n < 1, none of which are physically meaningful. Therefore, 
when n < 1, the solutions to the model equations break down physically in a finite time, and 
in that sense there does not exist a stable, moving reaction front. Nevertheless, it should be 
noted that Hilhorst et al. |jl5|, 16] have shown that well-defined solutions with free boundaries 
(at the diffusive scale x oc can exist when < ra < 1 . 



3 Derivation of the Asymptotic Similarity Solution 

3.1 Scaling of the Reaction Front 

The initial-boundary- value problem ([T0|) possesses no natural length or time scale, i.e. it 



is invariant under power-law "stretching transformations" |2^, and consequently in the limit 
t — » 00 the system is expected to approach an asymptotic similarity solution in which distance 
and time are coupled by power-law scalings Since reactant A diffuses while reactant 

B does not, the (presumably unique) point of maximal reaction rate r{a, h) = a^lf" moves 
in the —x direction toward the reservoir of reactant B. Therefore, an asymptotic similarity 
solution, if one exists, must involve a moving frame of reference centered on some point Xf{t) 
identifying the position of the reaction front at or near the point of maximal reaction rate 
(with dxf/dt < 0). Let x/(t) = — 2z/t°", where z/(g) > is a constant (akin to the "speed" 
of the front) to be determined self-consistently during the analysis, and consider an arbitrary 
coordinate stretching transformation in the moving reference frame, 

where w{t) = t° is the width of the reaction front indicated in Fig. |l|. (The factor of two is 
included only for algebraic convenience.) 

In the neighborhood of Xf{t), we also allow the magnitude of a{x, t) to vary with a power- 
law scaling, 

A{7],t) = f<a{x,t). (17) 

If 7 7^ 0, then another similarity solution far away from the reaction front (in the "diffusion 
layer" shown in Fig. |I| and defined below) will be needed to satisfy the boundary condition 
a(oo, t) = 1. This possibility that two regions with different asymptotically self-similar dynam- 
ics for a{x, t) could arise is suggested by the fact that there are two driving terms, representing 
diffusion and reaction, on the right-hand side of (|10a|) with different behaviors under stretch- 
ing transformations. On the other hand, there is only the reaction term on the right-hand 
side of (10b), so b{x,t) can exhibit only one type of asymptotic scale invariance. This is the 
main mathematical consequence of the physical fact that reactant B does not diffuse. Since 
6(— 00, 1) = 1, we consider the transformation 

B{r^,t) = h{x,t). (18) 

Note that the reaction term r{a,b) = a^h^ has the scaling, r = t'^A^B"", where f3 = rwy in 
the notation of Galfi and Racz iTH. 
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These transformations leave the governing equations in the form: 



t 



(m— 1)7 



dA 
'dt 



tim-lh-l (^^ + ^^^^ 

or] 

BE BE BE 

dt Br] Br] 

We now look for asymptotically invariant solutions 

BA 



^(m— 1)7— 2q! 



B^A 

Br]"^ 



-qA'^E^. 



(19a) 
(19b) 



A{r],t)^ A{r]), 
E{r],t)^E{r]), 



Br] 
BE 
Br] 



(r/,t)^^'(r/), 
M^E\i^) 



B^A 



{r^,t)~.A"{v) 



(20a) 
(20b) 



as t — >■ 00 with |?7| < 00 fixed. For consistency with our definition of the reaction front, we 
require that there is in each equation a dominant balance between the reaction term A^E"" 
and at least one other non-vanishing term. In order for time invariance to be attained, we 
assume that time- dependent terms in the transformed coordinates are negligible compared to 
the reaction term, i.e. 



liin t(— 1)7^ = 

t^oo Qt i^oo Bt 



0, 



(21) 



which is a precise statement of the assumption of "quasi-stationarity" |14 



A dominant balance with the reaction term in (|19a]) implies that at least one of the 
following three cases must be true: 



Case Al: 


(m — 


1)7 


-2a = 0, 


(m — 


1)7- 


1 


- a + cr < 0, 


(m — 


1)7- 


1 < 0, 


Case A2: 


(m — 


1)7 


-2a < 0, 


(m — 


1)7- 


1 


— a + 0" = 0, 


(m — 


1)7- 


1 < 0, 


Case A3: 


(m — 


1)7 


-2a < 0, 


(m — 


1)7- 


1 


- a + a < 0, 


(m — 


1)7- 


1 = 0. 



Likewise a dominant balance in (19b) requires that one of the following two cases must hold: 
Case Bl: 7717 — 1 = 0, m7 — 1 — a + cr<0. 
Case B2: 7717 — 1 < 0, m7 — 1 — a + cr = 0. 

There are only two combinations of these cases that are logically consistent: 
Traveling Wave Case: (A2, B2) a > 0, 7 = 0, cr = 1 + a, and 

Diffusing Front Case: (Al, B2) a = (m - 1)7/2, a = 1 - (m + 1)7/2, < 7 < 1/m. 

In the first case, we have a > 1, which implies that the reaction front advances at least linearly, 
e.g. as a traveling wave Xf ^ t, but in the second case, the front advances sublinearly, e.g. as 
a diffusing front x/ ~ t^/^. In both cases, note that the reaction order n of the static species B 
plays no role in the scaling behavior. The same conclusion is also true of the reaction order m 
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of the diffusing species A in the Travehng Wave Case, but m does affect the scahng exponents 
in the Diffusing Front Case. 

Consider the possibihty 7 = 0, which is only consistent with the Travehng Wave Case. In 
this case, a single asymptotic scale invariance is attained everywhere, and the equations for 
A{ri) and B{r]) are 

2auA' = 5o,,oA" - A'^B" (22a) 
2auB' = -qA'^B'' (22b) 

where 5x,y is the Kronecker delta. By combining these equations and integrating once using 
the boundary conditions behind the front, i.e. ^(oo) = l, ^'(oo) = and B{oo) = 0, we 
obtain 

2aiy{B + q) = q{2auA-6a,oA'). (23) 

Applying the boundary conditions ahead of the front, ^(— oo) = ^'(— oo) = and B{—oo) = 
1, to this equation then implies o"z/(l + g) = 0, which is a contradiction since a > and u > 
are needed for the reaction front to move at all (and g > 0). 

In this way, we are forced to consider at least two regions with different scale invariance 
if there is to be any hope of an asymptotic similarity solution. Since the second type of scale 
invariance is associated with the dominance of the diffusion term versus the reaction term in 
( |10aj ), it must occur only on the back (+x) side of the reaction front due to the reservoir of 



reactant A at infinity, a(oo, t) = 1 (see Fig. |l]). To describe the scale invariance of the diffusion 
layer, we postulate another power law W(t) = t^ for the asymptotic width of the diffusion 
layer. 

3.2 Scaling of the Diffusion Layer 

Since 6 a, there are two possibilities, each involving a singular perturbation w/W = t"~^: 
Infinitely Thin Reaction Front Case: 5 > a, w = o(H^), 
Infinitely Thin Diffusion Layer Case: 5 < a, W = o{w). 

Since chemical reactions are typically much faster than diffusion, the former case seems more 
reasonable on physical grounds, but we do not rule out the latter case a priori. In the Infinitely 
Thin Reaction Front Case, the reaction front is defined hj x — xj = 0{w) and the diffusion 
layer by W = 0{x — Xf), x > Xf, whereas in the Infinitely Thin Diffusion Layer Case, the 
reaction front is defined hj w = 0{x — Xf), x < xj and the diffusion layer hj x — Xf = 0{W). 
In both cases, we view the reaction front as representing the "inner problem" (with similarity 
variable \r]\ < oo) and the diffusion layer as representing the "outer problem" (with similarity 
variable C > defined below). The two regions are connected by asymptotic matching of the 



limits ?7 — > oo and C — ^ 0^ (described in the next section) ^ . 

In order to treat the outer problem, we transform the original equations using a new 
reduced coordinate with power-law scalings, 

X + 2i't'^ 

C= , AiC,t)=aix,t), BiC,t)=bix,t). (24) 
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(Another factor of two is included in r] = 2t^'^°'( for algebraic convenience. Note the use of A 
and B for the diffusion layer versus A and B for the reaction front.) Under this transformation, 
the equations take the form, 



t^^— -t^^-i^C^ + t^+^Vz.— = -^^-t^'A-B-, (25a) 
t^s^_^2S-i^^^_^^5+a-i^^^ = (25b) 

Seeking an asymptotic similarity solution, we assume that invariance is achieved in the 
transformed equations 

^(C,t)-A(C), |^(C,t)-A'(C), ^(C,t)-^"(C) (26a) 

%,t)-i?(C), ^(C,t)-i?'(C) (26b) 

assuming time-variations in ( p6a| ) are small relative to the diffusion term 

dA 

^hmt^^ — (C,t)=0, forC>0. (27) 

In order to obtain a different scaling from the reaction front, the reaction term must also not 
enter into the dominant balance 

limt2'5i((,t)™^((,t)" = forC>0, (28) 

t — ^oo 

a condition that we will check a 'posteriori for consistency in section ^. Note that this limit 
vanishes trivially for C, > v since we have already noted that 6(x, t) = for all x > 0. From 
(25b), this condition would imply dBjdC, = 0, which together with the boundary condition 
B{oo) = would imply B{() = 0. With the reaction term gone, one of the terms on the 
left-hand side of ( P5a| ) must balance the d'^A/dC,'^ term on the right side; if not, we would 
have d'^AjdC^ = 0, which cannot satisfy all of the boundary conditions. There are only two 
possible dominant balances: 

Case Dl: 5 + a - 1 = and 25 - 1 < 0, 

Case D2: 5 + a - 1 < and 25 - 1 = 0. 

The former case implies cr > 1/2 and hence contains the Traveling Wave Case (and not the 
Diffusing Front Case). With the scaling relations of case Dl, Eq. ( |25a) ) has the asymptotic 
form 

auA' = -A". (29) 
4 

The solutions to this equation exhibit exponential growth a.s ( —>■ oo, which is incompatible 
with the boundary condition A{oo) = 1. Therefore, we conclude that Case Dl, and hence also 
the Traveling Wave Case, is not consistent with the boundary conditions. At this point, we 
are left with Case D2 together with the Diffusing Front Case (Al,B2), which imply W ~ 
thus justifying the term "diffusion layer" for the region C > 0. 
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3.3 Asymptotic Matching of the Reaction Front and Diffusion Layer 

One more condition is needed to uniquely determine the scaling exponents, and it comes from 
asymptotic matching: The "outer limit" rj ^ oo of the inner approximation must match 
with the "inner limit" C 0^ of the outer approximation (because both are asymptotic 
representations of the same function). Unlike the more familiar case of boundary layers of 
ordinary differential equations however, our system of partial differential equations 

will require extra care for matching because the limit t ^ oo (with either ( or r] fixed) must 
be taken before the inner and outer limits. 

Since B{() = for all C > 0) only matching condition for b{x,t) is trivial, B{oo) = 0, 
but the matching conditions for a{x,t) are more subtle. Since 7 > 0, the concentration of 
species A approaches in the reaction front: a{x,t) = 0{t~^) as t — 00 with \ri\ < 00 fixed. 
Therefore, a boundary condition on the outer problem is A{() = 0, but unfortunately this does 
not provide a boundary condition on the inner problem. Instead, we must consider matching 
at the next (linear) order of Taylor expansion in the intermediate region: 

^ _ f Iff as t ^ 00 with < C < 00 fixed 

as t ^ 00 with Ir^l < 00 fixed " ^ ' 

Now requiring that the two intermediate limits match yields the final scaling relation: 

a + 7 = 5 (31) 

as well as the missing boundary condition on the inner problem: 

^'(00) = A'{0)/2. (32) 

(Note that ^(C) is already fully determined by matching at zeroth order.) This scaling relation 
(|3TD can be understood physically as expressing conservation of mass between the diffusion 
layer and reaction front 0]. Similarly, the matching condition (|3^ ) simply means that the 
diffusive flux entering the reaction front equals the flux leaving the diffusion layer. 

By examining all possible similarity solutions with power-law couplings of distance and 
time, we flnally arrive at a unique set of scaling exponents from cases Al, B2 and D2 and 

m-l m 1 1 

2(m + 1) m + 1 m + 1 2 

Therefore, after long times the reaction front itself "diffuses" according to x/(t) = — 2z/-\/t, 
where /^(g)^ is now interpreted as an effective diffusion constant for the front. Although for m = 
1 the reaction zone settles down to a constant width (a = 0), for m > 1 the front width grows 
in time (a > 0). In all cases the reaction front is "inflnitely thin" compared to the diffusion 
layer [a < 6). Note that as m increases, 7 tends to zero, meaning that the concentration 
in the reaction front does not decrease as quickly for higher-order reactions as it does for 
flrst-order reactions. The exponents a = and 7 = 1/2 for the case m = 1 were flrst obtained 
by Jiang and Ebner ^ based on physical arguments supported by Monte Carlo simulations 
and later discussed in an analytical context by Koza |l^, but to our knowledge prior to this 



work neither have the general expressions for m 7^ 1 been given nor have the scaling exponents 
been proven to be unique. With the scaling exponents and matching boundary conditions now 
determined, we proceed to solve the inner and outer boundary-value problems in the following 
sections. 
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3.4 Concentration Profiles in the Diffusion Layer 

In the diffusion layer from (25)-(^) we have 

- 2(C - iy)A' = A", A{0) = 0, A{oo) = 1. (34) 

The exact solution to this boundary value problem can be expressed in terms of error func- 
tions ||30l 

and is depicted in Fig. |^. Note that the dimensionless flux entering the reaction front 

■^'(^' - ^ - V?(l+erfM) 

is needed for asymptotic matching in (0). 

The effect of varying q = n' p°^/ {m' p°^) is easily understood in terms of the mathematical 
model. As q is decreased, the reaction front slows down since reactions in the front region 
remove species A much faster than diffusion can replenish it. In the limit g ^ 0, i.e. 0, 
the front comes to a complete stop, z/(0) = 0. For very small, but finite g' > 0, the concentration 
of diffusing reactant approximately obeys 

da 3'^ (X 

= a(x,0)=H(x), a(0,t) = 0, a(oo,t) = l, x>0 (37) 

ot ox^ 

(at least for t ^ z/~^ since the reaction front is stationary only for short times). This classical 
diffusion problem has the exact similarity solution 



a{x, t) = erf 



X 



{31 



which is precisely the z/ = curve in Fig. ^. From (^5[) note that even when the front has 
moved significantly {t 2> z/~^) the concentration still has the same shape, A{() ^ erf(C), in 
the (very slowly) moving reference frame, as long as z/(g) -C 1. 

On the other hand, as q is increased, the reaction term becomes progressively less important 
compared to the diffusion term in ( |10a| ). In the limit q — * oo, i.e. — * 0, we recover another 
classical diffusion problem (after sufficiently long times t ^ z/^^) 



da 

dt 



dx 



2 ' 



a(x, 0) = H{x) 



'00,t) = 0, a(oo,t) = 1, 



which has the exact similarity solution 



a{x,t) = - 



erf 



X 



2v^- 



(39) 



(40) 



Note that in the limit q ^ oo the reaction front instantly speeds off to — oo having consumed 
only a negligible amount of reactant A, resulting in a pure diffusion problem for a{x, t). Indeed, 
we will see below that z/(cx)) = oo. It remains, of course, to relate z/ and q. 

In Fig. 0(b), we see how the true asymptotic similarity solution A{() interpolates between 
the limiting forms (|38D and (0) as u goes from to oo, respectively. It turns out that for u > 2 
the asymptotic behavior of the original reaction-diffusion system is almost indistinguishable 
from (13) for t > 1 /4. 



12 



3.5 Diffusion Constant of the Reaction Front 



In the reaction front, we have a third-order system of nonhnear ordinary differential equations 

= ^"-^'"S" (41a) 

vB' = -qA'^B'', (41b) 

with four boundary conditions 

^(-oo) = 0, B{-oo) = 1, B{oo) = 0, and ^'(oo) = Ai, (42) 



where Ai{u) is known via (|36|) . Ahhough this boundary- value problem appears to be overde- 
termined, the fourth boundary condition is actually necessary to determine the unknown 
diffusion constant of the reaction front i^iq)- 

By comparing (19) and (41), some physical insight into the dynamics of the reaction front is 
gained. The concentration of diffusing reactant A is determined by a local balance of reactions 
and "steady state" diffusion and the concentration of static reactant B by a local balance of 
reactions and fictitious advection due to the translating reference frame. The latter balance 
reflects the special character of = 0: Since reactant B cannot diffuse to the front, instead 
the front must diffuse to it. This is no longer true if Db > (no matter how small), which 
explains why different scaling exponents arise in that case JI], These physical properties 
are manifested in the mathematical model by the fact that since it multiplies the highest 
derivative in the equations > is a singular perturbation. 

One integration of (41) is easy to perform and fortunately suffices to derive an exact 
expression for z/(g). Substituting (41b) into (f41aj) , integrating and applying the boundary 
conditions at = oo, we obtain 

iyB = q{Ai-A'). (43) 

Likewise enforcing the boundary conditions at r] = — oo, we find z/ = qAi. Substituting Ai 
from (^), we have 

uiq)=F-\q), (44) 

where 

F{x) = v^a;e^' [1 + erf(x)] , (45) 

which was first derived by Koza The function i>{q) is plotted in Fig. ^. 

The transcendental function F{x) cannot be inverted analytically, but limiting formulae 
can be derived. The Maclaurin series of F(x) is 



F{x) = V^x + 2x' + V^x' + \x^ + ^a;^ + ^x'' + ^x' + . . . , (46) 

3 2 15 

which can be inverted term by term to generate the Maclaurin series of 2^(5'), valid for small 

'^(^) = - + ^1 - • • • • (47) 

For large g, approximations such as, 

v{q) ~ 

can be generated by iteration. 



^ log (^-^^ - log2 - i log log ( ) , (4^ 
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3.6 Existence and Uniqueness of the React ion- Front Scaling Func- 
tions 



With the results of the previous section, the inner boundary-value problem is reduced to a 
nonlinear, second-order equation for A{ri): 

A" = A'^il - A'/AiT, Ai-oo) = 0, ^'(oo) = A^. (49) 



Once this system is solved, B{ri) is recovered from B{r]) = 1 — A'{r])/Ai. Note that ( 149|) 
is invariant under translation t] t] — rjo, where the arbitrary constant rjo sets the precise 
location of the reaction front. Since t]o depends on the exact initial conditions, however, it 
cannot be determined by considering only the long-time asymptotic limit as we have done 
here. 

Since the second-order equation ( ^91) is autonomous (i.e. r] does not appear), it is useful 
to consider the "Lie diagram" or "phase plane" [^, |29| of trajectories in the {A, B) plane 



parameterized by r], as shown in Fig. ^. By studying properties of the phase plane, it is 
straightforward to prove the existence and uniqueness of solutions if and only if m, n > 1 and 
^1 > 0. With the change of variables 

s = 4"^-^)/('"+^)r/, u{s) = B{r]), v{s) = Af^^'^^'Uiv), (50) 

we begin by transforming ( |49| ) into a system of first-order equations 

u' = -v'^u'' (51a) 
v' = l-u (51b) 

with boundary conditions f (— oo) = and u{oo) = 0, or equivalently u{—oo) = 1 and 
v'{oo) = 1. 

There is a unique fixed point at {u,v) = (1,0) corresponding to the region ahead of the 
reaction front which contains only the static reactant B. This is the starting point (s = — oo) 
of any trajectories that satisfy the boundary condition v{—oo) = 0, so our task is to identify 
and follow any unstable manifolds leaving (1,0) to see if they satisfy the other boundary 
condition u{oo) = 0. If m = 1, then the equations can be linearized about the fixed point, 
and (1,0) is a hyperbolic saddle point with an unstable manifold in the (1,-1) direction and 
a stable manifold in the (1, 1) direction, as shown in Fig. ^(a). 

If m 7^ 1, then the stable and unstable manifolds are degenerate at linear order and form 
a cusp oriented in the (0, 1) direction, as shown in Figs. §(b) and (c). The nonlinear stability 
of the fixed point can be determined by noting that f" ~ as -u 1. This equation has 
solutions satisfying the boundary condition v{—oo) = if and only if m > 1. A stable reaction 
front does not exist if m < 1 because the concentration of diffusing reactant A would become 
negative ahead of the front. Therefore, since m > 1 implies a = (m — l)/2(m + 1) > 0, the 
front width w{t) ~ either stays the same (for m = 1) or increases (for m > 1) but cannot 
decrease in time. For m > 1, we integrate v" ~ v'"^ once and substitute into (51b) to obtain 
the separatrices in the upper half-plane (f > 0): 



« ~ 1 ± as {u, v) ^ (1, 0+), (52) 
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where the upper sign corresponds to the stable manifold [r] oo) and the lower sign to the 
unstable manifold (77 — —00). 

Let us briefly consider trajectories in the lower half-plane [v < 0) in the neighborhood 
of the fixed point (m ~ 1). Of course, such trajectories are not physically allowed, but it is 
satisfying to prove that the model equations exclude such possibilities. If m is either irrational 
or a rational number of the irreducible form fci//c2 where k2 is even, then such trajectories do 
not exist because in that case (and hence u') would not be a real number. If m = ki/k2 
where ki is even and ^2 is odd, then the direction field {u', v') is an even function of f , which 
in light of (^) implies that there are no other separatrices in the lower half-plane and that 
trajectories merely circle the fixed point, as shown in Fig. ^(b). Finally, if m = ki/k2 where 
both ki and /c2 are odd, then u' is an odd function of v (while v' is even), and (1, 0) is a saddle 
point. In this case, Eq. (^) also describes separatrices in the lower half-plane with the lower 
sign corresponding to the stable manifold and the upper to the unstable manifold, as shown in 
Fig. |^(c). As it leaves the fixed point, this branch of the unstable manifold enters the region 
(m > 1, < 0), throughout which m' > and f ' < 0, and thus it heads off to -u = = — 00 and 
cannot satisfy the other boundary condition, as shown in Figs. §(a) and (c). Therefore, any 
solutions must lie entirely in the first quadrant of the phase plane (m > 0,f > 0). 

In this way we are left with only one possible solution, which leaves the fixed point along 
the unstable manifold of (|52| ) and enters the region defined by < m < 1 and f > 0, throughout 
which m' < and v' > 0. Since the f-axis {u = 0) is itself a trajectory, which cannot be crossed, 
this candidate solution must reach an asymptote u{oo) = Uo, for some constant < Mo < 1. 
However, it is clear from ( |51a|) that Mq = is the only possible asymptote, which implies that 



the trajectory (if it exists) must satisfy the other boundary condition u{oo) = 0. To check the 
existence of this solution in the limit s 00, note that f ' ~ 1 which implies u' ~ — (s — So)™^" 
for some constant So with m > 1. Solutions to this equation satisfying u{oo) = exist if and 
only if n > 1. A stable reaction front does not exist if n < 1 because the concentration of 
static reactant B would be negative behind of the front. 

3.7 Concentration Profiles in the Reaction Front 

Although solutions to (^) exist for m,n > 1, they are not easily expressed in terms of 
elementary functions. The exact trajectories in the phase plane, however, can be obtained. 
The ratio of (51b) and (|51a]) yields a separable, first-order equation for v{u) 



? = ^ (53) 

which can be integrated to obtain the one-parameter family of trajectories (for u > 0) plotted 
in Fig. I 

{ci + u — \ogu if n = 1 

C2 + \ogu + u~^ ifn = 2 (54) 

Cn - (n - 2)^2-" + {n- 1)-^^^-" if n > 1, n 7^ 2 

indexed by the real number c^. 

Applying the boundary conditions v{u = 1) = (which determines c„) and v{0) = 00 
(which selects the positive branch when v{u) is multivalued), we arrive at the exact phase- 
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plane trajectories (in the region < m < 1, f > 0) of the solution to the inner problem: 



[(m + 1) (M-l-logM)]'/^'"+'^ 
[(m + 1) (?i-i-l + logu)]'/^™+'^ 

(m+l)(l-(n-l)«2-n + („_2)«l-")"| 
(n-l)(n-2) 



if n 
if n 



1 
2 



(55) 



if n > l,n 7^ 2 



which is an algebraic equation v = gm,n{u) relating Airj) and B{vi) via (|50D. This equation is 
transcendental, but in some cases it is easily solved for u = g^\{v), e.g. for n = 3 we have 



u 



2(m + l)-if;'"+i - 1 
2(m + 1)-^-^+^ - 1 



(56) 



Note that fi'm^nlO) = 1 and 5'm^„(cx)) = 0. 

By substituting (|55D into (51b) we arrive at a first-order equation for virf) 



(57) 



without any boundary conditions (because the conditions v{—oo) = and v'{oo) = 1 are 



automatically satisfied). Since 
expressed in the form: 

A{v) = 
B{v) = 

where 



is separable, the solution to the inner problem can be 



■^1 ""m.n {•^1 



,(m-l)/(m+l) 



iv - Vo)) 



9m,n 



hr 



ds 



9v 



(58a) 
(58b) 

(59) 



The precise location of the reaction front is set by choosing v{so) = Vo > for some constant 

;ed to the solution of two 
). In practice, however, it 



So = Ai^ ^^^^"^^^''rjo. With these results, the inner problem is reduced to the solution of two 



algebraic equations for g^\{v) and h^^{ri) and one quadrature 
is simpler to integrate (pQf ) directly. 

Numerical solutions of the rescaled inner problem 



v[ — oo 



0, 



V [OO] 



1 



(60) 



are obtained by a shooting method, with the results shown in Fig. ^(a) for m = n = 1 and m = 
n = 2. The position of the front is chosen such that f (5) = 5. The static reactant concentration 
u{s) and the reaction rate density v'^u^ are shown in Figs, ^(b) and (c), respectively. Note 
that the concentration fields decay to their asymptotic values as \ri\ oo more slowly as m 
and n are increased above unity, a phenomenon that we explore analytically in next section. 



3.8 Localization of the Reaction Front 

The width of the reaction front varies in time according to w(t) ~ where a = (m— l)/2(m + 
1). On the scale of the diffusion layer width W{t) ~ t^^^, the reaction front is "localized" 
after long times because a < 1/2. Note that the overall localization of the front w{t)/W{t) 
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is controlled by m (the reaction order of the diffusing species A), but we now show that both 
m and n (the reaction order of the static species B) affect localization on the scale of w{t). 
Specifically, we derive the spatial decay of the scaled reaction rate 7l{ri) = A{i])"^B{r])^ in 
terms of the inner similarity variable r] —>■ ±00 (see Fig. ^(c)). The actual reaction rate decays 
uniformly to zero in time, r ~ t^^TZ^rj) with j3 = m/ {m + 1), but here we are only concerned 
with the shape of 71(1]). 

Ahead of the reaction front in the limit rj —00, we have A!' ~ from (^), which can 
be integrated to obtain the decay of concentration fields: 

A( \ I ^2e-l''l ifm = l 

l-^(^) - I ^3|r^|-(-+i)/(-i) if,„>l (61b) 

f ^'"e~™l''l if m = 1 

^(^) ^ I ^-|r/|-2W(r.-i) ifm>l (6^^) 

where A2 and ^3 are constants. Note that the localization of A{ri) and TZ{ri) ahead of the 
front is entirely controlled by the reaction order m of the depleted reactant (which is the 
diffusing species A). There is a transition from an exponential decay for m = 1 to a slower 
power- law decay for m > 1. 

Next we consider localization of TZ{r]) behind the reaction front in the limit 77 — 00. From 
asymptotic matching with the diffusion layer we have already derived A{r]) ~ Airj. The 
asymptotic decay of B{ri) and 7^(77) is obtained by integrating AiB' ~ —{Airi)"^B"': 



^ ^ 1 i3ir7-("^+^)/("-i) if n > 1 ^ ^ 

^ r ATB^v"" exp [-AT-'v'^+'/im + 1)] if n = 1 
^ \ ^™i3«r/-("^+")/("-i) if n > 1 ^ ' 

where Bi is a constant. Once again, a higher reaction order n for the depleted species (which 
is the static species B) broadens the front: There is another transition from exponential decay 
for n = 1 to a power-law decay for n > 1. Note, however, that increasing the reaction order 
m of the diffusing species A contracts the back side of the front. 

The fact that TZ{ri) has a fairly broad, power-law decay for n > 1 toward the diffusion 
layer should cause concern since we have previously assumed in (pSf) that the reaction term 
is negligible in the diffusion layer. Fortunately, however, for all m,n > 1 the decay of B{ri) is 
just fast enough to satisfy ( pSf ) in the intermediate region where B{0 < C ^ 1) ~ B{ri ^ 1). 
From ( |62a| ) along with A{C,t) ~ A{C), B{ri,t) ~ B{ri), rj = 2(tV{'n+i) and 5 = 1/2, we have 

t^^A(C,trB(C,tr = { \ ,L n u\ 63 

^ V^' ^ I Q f^-n{m+l)/(n-l)^-n/in-l)\ if n > 1 

as t ^ oo with < C ^ 1 fixed, which verifies (^Sj) in the intermediate region for any m,n > 1. 
However, we now prove that ( pSf ) actually holds throughout the diffusion layer for all ( > 0. 
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4 Transient Decay in the Diffusion Layer 



4.1 Assumption of Quasi-Stationarity 

The analysis of long-time asymptotics in the previous section rests on two basic assumptions: 
(i) "scale separation" , given by 

< Mm A{7],t)"'B{C,t)'' < oo for|77|<oo (64a) 
limt-i(C,t)'"5(C,t)" = forC>0, (64b) 

t — ^oo 

and (a) "quasi-stationarity" , given by (20) and (26) along with 

dA 

hmt(™~i)/("+^) — (r/,t) = for|r7|<oo (65a) 

^limt"^/("+i)^(r/,t) = for|r7|<oo (65b) 
dA 

lim t ■ — (C,t) = for C > 0. (65c) 

t^oo (jt 

Assumption (z) states that two spatial regions with disparate length scales, the reaction front 
and the diffusion layer, arise where the reaction term on the right-hand side of the governing 
partial integro-differential equation ( |T^ ) or (|T5|) is, respectively, either comparable to or dom- 
inated by the diffusion term. Assumption (ii) states that, when viewed on scales appropriate 



for each region, the solution to the initial-boundary- value problem (10) approaches an asymp- 
totically self-similar form, which is suggested by the fact that there is no natural length scale 
in the problem. 



These ubiquitous assumptions |T], |T3|, Q have been rigorously justified ||TT], in the 



special case of a perfectly symmetric {Da = Db, p\ = Ps)) thus stationary {y = 0), 
reaction front involving two diffusing reactants with certain kinetic orders {m = n = 1 and 
m = n > 3). To our knowledge, a similar mathematical validation of these assumptions has 
not been given for the general situation of a moving reaction front with arbitrary kinetic orders 
for either one or two diffusing reactants. For one static reactant, however, convergence at the 



diffusive scale has been rigorously established by Hilhorst et al. |TH], even with a very general 



reaction term |16] 



In this section, we prove the more modest result that quasi-stationarity implies scale 
separation, i.e. ( |65D implies (64). (It suffices to show (64b) since (64) follows from the 
definition of rj in section |0|.) Although this analysis justifies a posteriori the assumption 
of scale separation in our fairly general situation {Db = 0, z/ > 0, m,n > 1), it more 
importantly reveals the transient decay to the asymptotic similarity solution in the diffusion 
layer. Specifically, we derive exact formulae for the asymptotic decay of the reaction-rate 
density and static-reactant concentration in the diffusion layer. 



We begin by precisely stating our assumptions related to (pq). From ( |65a| ) and its 
consequence ( |61a| ), we conclude that the diffusing reactant concentration vanishes on the 
scale W{t) = \/t ahead of the front (where there is no diffusion layer) since A{(,t) ~ 
A{2(t^^^"^^^^) ^ as t ^ oo with C < fixed. This result can be combined with ( P6a|) 
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to obtain a statement of quasi-stationarity on the scale W{t): 

BA 

A{C,t)^Ai()Hi() and _((, t) A'(C)i/(C) (66) 
as t ^ oo with C 7^ fixed. We have aheady derived the exact form of the similarity function 



A(Q in (|35D as a consequence of neglecting reactions (64b). To avoid a circular argument, 
however, we must now establish (64b) without using (pSf), thus giving a postenon justification 
for the latter equation. Throughout section ||our only assumptions about A{() are ^4(0) = 
and A'{0) > 0. These properties follow from matching with the reaction front, where a — 



[i.e. 7 > 0) follows from quasi-stationarity, as shown in section 3.1 



4.2 Direction of the Diffusing- React ant Flux 

Let us prove that A{() is strictly increasing in the diffusion layer, A'{() > for all C > 0, as 
a consequence of (64). Combining ( |65c| ) with ( |25a| ), we have 



OA A 

2(C-'^)^ ~ ^-i^"^" as t ^ oo with C > fixed (67) 



which is easily integrated once using an integrating factor. 



^(C,t)-e-(^-^)^ 



e^' — (0,t)+t / i(e,t)'"5(e,t)'^e«-^)'t/e 



(68) 



Since the second term on the right-hand side is non- negative for alH > 0, we can pass to the 
limit t — > cxD for any fixed C > to obtain the desired bound 

A'iC) > A' {0)6"'-^'^"'^' > 0, (69) 

which expresses the physical fact that everywhere in the diffusion layer a nonzero fiux of the 
diffusing species is directed toward the reaction front (at sufficiently large times). 

4.3 Decay of the Static- React ant Concentration 

We now prove that b{x, t) vanishes asymptotically in the diffusion layer as a consequence of 
quasi-stationarity, which implies 

5(C, t) B{C) = H{-C) as t ^ oo with C 7^ fixed. (70) 

For C < 0, this follows from (61b) since there is no diffusion layer ahead of the front, and 
therefore B{C,t) ~ 8(2(1^/^'''+^^) 1 = B{C) as t ^ oo with C < fixed. Likewise, in section 
2?1| , we have already established (^) for C > since B{(,t) vanishes there identically for all 



times. Therefore, it only remains to prove that -B(C) = for < ( < u. 

In light of the expression for b{x,t) in (p!2|), the definition of (j)rn{x,t) in (p!3| ) and the 
restriction n > 1, it suffices to show that ^m{(,t) = (f)rn{x,t) — >cx)ast^ooforO<C<i^ 
fixed. Using ( = u + x/2v^ (since 5 = cr = |), we transform (p^ix^t) into the diffusion-layer 
coordinates (x, t) ^ iC^t) 

*m(C, t) = /^(^ (2v^(C - r^dr, (71) 
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and express this in terms of the diffusion-layer scahng function a(x, r) = yl(z/ + x/2y^, t) 



<i>U(,t)= f'AlJ-{(-u) + u,T] dr. 



(72) 



It is convenient to work with the partial time-derivative of $m(C;^) given by the Leibniz rule: 

Focusing on the region < ( < u, we make the transformation ^ = \/t/T{C ~ ^) ^ ^i 

di, (74) 



dt ' J-oo dC \ V'^ ~ 

and pass the limit t ^ oo inside the integral to obtain 



iim — — 

t-*co dt 



(75) 



where the lower limit of integration follows from ( |66D since A{() = A'{() = for < 0. 
This step is justified by the Dominated Convergence Theorem because, by virtue of (|66D , 
there exist constants M,to > such that the integrand in ([7^ ) is bounded for all t > to by 
M/((^ — i/)^, which is integrable on {—oo, (), if ( < u. 

Since ^4(0) = is required by matching between the two regions of quasi-stationarity, 
Eq. (|75D can be written in the form d^m/dt ~ fm{C), where 



/-(C) 



C dA" 



dC 



■(0 



di. 



(76) 



Note that /^(C) > for < C < since A'{Q > in this region, as shown in (^). Therefore, 
with an integration of (p6|), we arrive at the desired result 

*m(C) ~ fmiQt as t ^ oo with Q<C<v fixed, (77) 

thus completing the proof that 5(C) = for ^ > 0. 

By substituting (fTzD into ([T^), we obtain the transient decay of B{(,t) in part of the 
diffusion layer where the reaction front has already passed {0 < ( < u): 



if n 



[g(n - l)/„(C)r'/^"~'^ ifn>l 



(7^ 



Note that B{(,t) vanishes with exponential decay if n = 1 and with a power- law decay if 
n > 1. Therefore, by measuring the asymptotic decay (either exponential or power- law) of 
the static-reactant concentration in the diffusion layer, the reaction order n could in principle 
be inferred from experimental data (although such measurements are difficult in practice [||). 
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4.4 Decay of the Reaction-Rate Density 



From (^) we easily obtain the asymptotic decay on the reaction rate density in the diffusion 
layer as t ^ oo with C > fixed: 



tA{c,trB{c,tr 



Aicr [qin - 1)/™(C)]" 



if n = 1 
if n > 1 



(79) 



' n-l f n-1 

which establishes (64b). The reaction term in the diffusion layer has previously been neglected 



based only on physical intuition [14|, but here we have given a mathematical justification. 



4.5 The Decay Time When n = 1 

Since the reaction term vanishes sufficiently fast in the diffusion layer to justify a posteriori 
the analysis in section the exact expression for A{() from (^) may be substituted into ([76|) 
to evaluate the function fm{C)- If ^ > 1, then fm{() affects the transient decay in (|78D and 
( |79D only as a multiplicative prefactor in a power law, which would be difficult to measure in a 
real experiment. If n = 1, however, then fm{C) sets the characteristic time Tm{C)~^ = ifmiC) 
of an exponential decay, which is easier to measure experimentally. Therefore, we now derive 
an exact expression for the decay time ti{Q (0 < C < z/) in the case m = 1: 



n(C) 



-1 



9/i(C) 

c 



1 - 



2{y - C) - C) - z/e*^'-^""^^' 
+v^z/e'^' [l + 2(z/ - C)'l ■ [erf(z/) - erf(i/ - C)] 



^0) 



Note that ti(0) = oo in the vicinity of the reaction front (C = 0) because (^) and (|78| ) no 
longer hold. Within the diffusion layer, the decay time is a decreasing function of distance C, 
away from the reaction front, as shown in Fig. ^ These results may be used to infer reaction 
orders and perhaps even kinetic constants in diffusion-limited corrosion experiments from 
transient decay measurements of the reaction-rate density in the diffusion layer H]. 



5 Uniformly Valid Asymptotic Approximations 

In the previous two sections we have argued for the existence of a unique asymptotic simi- 
larity solution (up to an unknown constant r^o) contingent upon certain "quasi-stationarity" 
conditions, which are likely to be satisfied for the specified initial conditions (see below). This 
solution, valid after long times, consists of two different asymptotic approximations for a(x, t), 
the concentration of the diffusing reactant A, which reflect the different couplings of length 
and time in the reaction front and the diffusion layer. A single asymptotic approximation for 
a(a;, t) that is uniformly valid across all space is obtained by adding the two contributions from 
the reaction front (the inner region) and the diffusion layer (the outer region) and subtracting 



the overlap (from the intermediate region) ||2^, |28 



a(x, t) ~ \A{r] - r]o) ~ Ai ■ {r] - r]o)H{r] - r]^)] t"^/('"+^) + A(C)i/(C), as t ^ oo for all x. (81 
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where the reaction-front and diffusion-layer similarity variables are 



V[x,t) = — -T— ^ — -rz (82a) 

' ^ +(m,-l)/2(m,+l) ^ ' 

and a^,t) = ^^^17^, (82b) 

z/(g)^ is the diffusion constant of the reaction front (see and Fig. Ai = i^{q)/q is 
a constant proportional to the diffusive flux entering the front, A{r]) is the reaction- front 
similarity function (see ( ^8a| ) and Fig. A{() is the diffusion-layer similarity function (see 
(|35D and Fig. 0) and rjo is an undetermined constant depending upon the initial conditions 
as well as the precise definition of the reaction-front location. The uniform approximation 
has been determined analytically up to the solution of two algebraic equations (58) and one 
quadrature (|59D . 

A subtle point in the construction of this uniformly valid approximation is that shifting 
the position of the front hj rj t] — rjo does not affect matching with the diffusion layer 
because in that case C '"^ (^7 ~ r]o)t~^^^"^'^^y2 ~ r]t~^^^"^~^^^ /2 = In other words, because 
the reaction front is "infinitely thin" compared to the diffusion layer, translating its similarity 
variable by a constant 7]o, or any other function of time that is o(t^/^), does not require that 
the diffusion-layer similarity variable ( be shifted as well. 

The situation for b{x,t), the concentration of the static reactant B, is much simpler. By 
comparing the asymptotic bound on b{x, t) in the diffusion layer given by ( [78| ) with the tail 
of the reaction-front approximation given by (|62a|) with r] = 2(t^^^'^~^^\ we see that that the 



asymptotic behavior of b{x, t) is identical in the two regions. Therefore, 

b{x, t) ~ B{ri — rjo), as t ^ oo for all x (83) 

is a uniformly valid approximation, where B{ri) is the reaction-front similarity function given 
by (58b). 

At this point the initial conditions have not yet entered the analysis except in (|l2l), which 
only influences the prefactors of the transient-decay formulae in section ^. Therefore, the 
asymptotic similarity solution is universal up to a constant shift of the reaction front by rjo for 
some broad set of initial conditions which presumably contains ( |10d| ). In general, this "uni- 
versality class" of initial conditions leading to the same asymptotic similarity solution (up to 
different values of rjo) is expected to be attained whenever the initial reaction-rate distribu- 
tion r(x, 0) = a(x, 0)™6(x, 0)"" is sufficiently well localized and the reactants are sufficiently 
well separated. This class surely contains all initial conditions for which r(x, 0) has compact 
support, e.g. r{x, 0) = for x 7^ in ( p.Od| ), or exponential decay, e.g. r(x, 0) < Me-l^l/^° for 
some M, Xo > 0, but perhaps not slower power-law decay. 



6 Discussion 

In this article we have studied the long-time asymptotics of solutions to the initial-boundary- 
value problem of (|^)-(6), which is a generic mean- field model for the corrosion of a porous 
solid by a diffusing chemical. We have derived a uniformly valid asymptotic approximation 
(|8TD-(^) consisting of matched similarity solutions in two distinct regions, the reaction front 



22 



and diffusion layer, each possessing different power-law scaling behavior. The existence and 
uniqueness of the similarity functions and the scaling exponents have been established if and 
only if m, > 1, and through an analysis of transients in the diffusion layer the asymptotic 
scale separation has been shown to follow from the assumption of quasi- st at ionarity. Since 
quasi-stationarity has been observed in recent experiments on the corrosion of ramified elec- 
trodeposits [Q, the present analysis therefore suffices to establish the theoretical predictions 
of the mean-field equations for at least one particular corrosion system. Although the case 
considered here {m,n > 1, g 7^ 1) is more complicated, it would be useful to perform a rig- 



orous transient analysis along the lines of Schenkel et al. |TT|] (who considered only the case 
Da = Db ^ 0, q = ^ and m = n = 1). Nevertheless, we have at least provided a firm 
mathematical justification for the scale separation between the diffusion layer and reaction 
front. 

In this work we have paid special attention to the effect of higher-order reactions (m, n > 1). 
First of all, the scaling exponents vary with the reaction order m of the diffusing reactant in 
precisely the same way as they do on the sum m + n in the case of two diffusing reactants 0, 
as shown in Table |l[ Moreover, the spatial localization of the reaction rate r{a,b) on each 
side of the front depends primarily on the reaction order of the depleted reactant: As the 
appropriate reaction order is increased from unity, the spatial dependence of the reaction rate 
away from the front changes from an exponential decay to a progressively broader power- 
law decay. Similarly, the temporal decay of the depleted (static) reactant concentration in 
the diffusion layer depends sensitively on its reaction order, undergoing a transition from 
exponential to power-law decay (in time) as n is increased from unity. These properties may 
have general relevance for more complicated multi-component reaction-diffusion systems. 

Other qualitiative features of our analysis that might have more general applicability are 
the dominant balances in the reaction front, where the concentration of a diffusing reactant is 
determined by a balance between reactions and "steady state" diffusion, i.e. a mobile species 
diffuses slowly to the front where it immediately reacts. On the other hand, the concentration 
of a static reactant is determined by a balance between reactions and fictitous advection 
due to the moving reference frame, i.e. since the static species cannot diffuse to the front, 
the front must diffuse to it. These guiding principles might help simplify more complicated 
reaction-diffusion equations for which asymptotically self-similar solutions do not exist. 



Epilogue 

After the writing of this article, a referee pointed out some interesting similarities between our 
analysis of chemical reaction fronts and various existing studies in combustion theory [^2|, |3^ . 



Indeed, the equations for combustion waves introduced by Zeldovich and Frank-Kamenetzki |^ 
which have since been generalized by many authors (e.g. Matkowsky and Sivashinsky ^^^), 
bear some resemblance to (2) since they describe the diffusion of a fuel substance coupled to 
the diffusion of heat. In combustion theory, however, the usual reaction term is quite different 
from d^) because it involves exponential Arrhenius temperature dependence, and the initial 
and boundary conditions also differ from those considered here. As a result, combustion waves 
tend to exhibit qualitatively different behavior from reaction-diffusion fronts. (For example, 
simple flame fronts have constant width and constant velocity.) Nevertheless, combustion 
waves exhibit multiple scales analogous to the reaction-front and diffusion layers described 
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here, which have also been analyzed using matched asymptotic expansions |^ (although not 
in the dynamical setting of this work). The idea of matching derivatives between the inner 
and outer regions actually appears to have its origin in the pioneering paper of Zeldovich 
and Frank-Kamenetzki [^], in which the velocity of a simple flame front is determined by 



analyzing a single-component equation like (§). In hindsight, it is somewhat surprising that 
the recent parallel literatures on two-species reaction-diffusion fronts and combustion waves 
have developed quite independently of each other, without any cross-references (at least, none 
to our knowledge). It is hoped, therefore, that this paper will initiate the "diffusion" of ideas 
between these two mature but related disciplines. 
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(i) Da>0, Db^O 
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2{m+l) 


m+1 
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m+n— 1 


m+n 
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2(m+n+l) 


m+n+1 


m+n+1 



Table 1: Comparison of the two possible sets of scaling exponents for asymptotic 
similarity solutions to the onc-dimcnsional, mean-field reaction-diffusion equations 
for two initially separated reactants. The width of the reaction front is given by 
w ~ i", the magnitude of the reaction rate by it! ~ and the concentration of a 
diffusing species in the reaction front hy a t~^. 
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Figure 1: Schematic diagrams (on the left) and concentration sketches (on the right) 
showing (a) the diffusing reactant A initially separated from the static reactant B and 
(b) the formation of a reaction front propagating into the region rich in species B at 
some later time leaving in its wake a mixture of species A and the reaction product 
C. Various quantities discussed in section |^ are also indicated, e.g. the position of 
the reaction front Xf{t) and the widths of the diffusion layer W{t) and reaction front 
wit). 
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Figure 2: The asymptotic similarity function in the diffusion layer, a{x,t) ~ 
A{C) where C = + x/2y/i shown for u = 0.0, 0.5, 1.0, 1.5, 2.0 versus C in (a) 
and versus C — in (b). The limiting shape CTi(x/2\/t) corresponds to u = 0. 
The other limiting shape [l + eTf{x/2\/t)]/2 as z/ — oo is plotted as the dashed 
line in (b), but it is almost indistinguishable from the i> = 2 curve. 
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Figure 3: The exact dependence of (the square root of) the dimensionless 
diffusion constant of the reaction front, on the parameter q = n' p\/ m' p% from 
Eq. (H. 
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Figure 4: Phase-plane trajectories for the inner (reaction- front) boundary- 
value problem from Eq. (|55D labeled by the constant Cn for m = n = 1 in (a), 
m = 2,n = 1 in (b) and m = 3, n = 1 in (c). In each case, the solid lines are 
separatrices emanating from the fixed point (1,0) of which the solution to the 
inner problem corresponds to the unique curve connecting (1,0) and (0,oo). 
Arrows indicate the direction of increasing s. 
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Figure 5: Structure of the reaction front obtained by numerical solutions of 
Eq. ( |60D with So chosen to set f (5) = 5. Profiles of (a) v{s) = A{ri)Ai 
(b) the static reactant concentration u{s) = B{ri) and (c) the reaction rate 
density u{s)"^v{s)^ = TZ{r))Ai'^'^^^'^^^^ versus s = are shown for 

m = n = 1 and m = n = 2. Arrows indicate that the reaction front propagates 
from right to left. 
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